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A perturbative method is developed for calculating the effects of recurrent synaptic interactions 
between neurons embedded in a network. A series expansion is constructed that converges for 
networks with noisy membrane potential and weak synaptic connectivity. The terms of the series can 
be interpreted as loops of interactions between neurons, so the technique is called a loop-expansion. 
A diagrammatic method is introduced that allows for construction of analytic expressions for the 
parameter dependencies of the spike probability function and correlation functions. An analytic 
expression is obtained to predict the effect of the surrounding network on a neuron during an 
intracellular current injection. The analytic results are compared with simulations to test the range 
of their validity and significant effects of the the recurrent connections in network are accurately 
predicted by the loop-expansion. 
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I. INTRODUCTION 

Recurrent connectivity is a ubiquitous feature of bio- 
logical neural networks, particularly in the central ner- 
vous system of vertebrates. This recurrent structure has 
made mathematical analysis difficult because of the non- 
linearities that arise when the action of a particular neu- 
ron affects itself through its action on neighboring neu- 
rons 0. In addition, interpretation of electrophysiologi- 
cal data can be misleading because the activity of neigh- 
boring neurons to the recording site are typically affected 
by the experimental stimulus @> 0> HI • 

Mathematical methods that can quantify the effects 
of neighboring neurons on the response of a single neu- 
ron embedded in a network would greatly improve the 
predictions of neural network models and ease the com- 
parison with system-level electrophysiological recordings. 
The development of theoretical methodology has led to 
the study of idealized systems such as syn-fire chains |^ , 
networks with correlated inputs [3,|8j, and large networks 
of identical neurons The goal of this project is to 

develop analytic methods to improve quantitative tests 
of electrophysiological hypotheses in order to set stricter 
limits on our theories of brain function. 

Quantitative tests require mathematical methods that 
separate different dynamic modes and can be refined to 
dynamics of the particular system of interest. We shall 
separate the dynamical modes of recurrence into coherent 
and asynchronous modes [12| . The coherent mode is 
characterized by strong synaptic connections that cause 
the activity of one neuron to dramatically affect the ac- 
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tivity of another. This dynamical mode is essential to 
central pattern generators. A simple example is the half- 
center oscillator, consisting of a pair of symmetrically 
coupled neurons that are strongly coupled by inhibitory 
synaptic currents. When one neuron fires, the other neu- 
ron is silent. If each neuron expresses intrinsic conduc- 
tances to terminate a burst, and to induce firing once 
they are released from inhibition, then the neuronal pair 
will oscillate. The coherent mode is quantified by the cor- 
relation function that is nearly ±1, for synaptic couplings 
that are either excitatory (+1) or inhibitory (-1). 

The present article investigates the effects of recurrence 
in the asynchronous state when the joint-correlation be- 
tween neurons is small but does not vanish. We quantify 
the parametric dependence of neural-response variables 
in the asynchronous state on the synaptic strength, the 
spike probability, and the presence of noise in the sys- 
tem. We also show that the influence of neurons on their 
neighbor can be considerable in the asynchronous state, 
suggesting that such corrections should be included in 
our predictions about the results of physiological exper- 
iments, even when no coherent state dynamics are obvi- 
ous. 

Another interesting result found in this analysis is that, 
under physiological conditions found in the vertebrate 
central nervous system, the effects of recurrent connec- 
tivity propagates across few synaptic contacts. Thus, if 
distance is measured by the number of synaptic junc- 
tions separating neurons, effects of distant neurons can 
be safely neglected in models of biological neural net- 
works. Thus, large-scale simulations of biological neural 
networks could be simplified to smaller networks that 
produce quantitatively similar results. 

In the following two subsections, we introduce our neu- 
ral modeling methods and establish our notation. The 
next section investigates the effects of recurrence on the 



2 



|Nois« i> 



Spike output/'iM ."'"'\t — »""'!'f — 14 




Background inpul ('((„.(,.' 



FIG. 1: Recurrent network model. Spike response neurons 
are coupled by recurrent synaptic connections with weight 
w. The spike probability, P(t n ), of each model neuron is 
a function of the synaptic connections, a background input, 
U(t n ), and internal noise, a. 



simple model of a pair of neurons. We show how to cal- 
culate the correction to the spike probability and mem- 
brane potential due to recurrent synaptic connections. 
The method considers recurrent connections as a pertur- 
bation to the background activity of the mean field result 
|13| . The terms of the perturbation series are weighted 
by powers of the synaptic strength multiplied by the in- 
verse of the noise in the system. Thus, the perturbation 
expansion converges quickly for systems that are noisy 
and weakly coupled. 

We compare the results of the perturbation expan- 
sion to numerical simulations to test the validity of the 
method and quantify the limits of this approach. Other 
experimentally measurable quantities are also calculated: 
the response to a stimulus, the auto-correlation function, 
and the joint-correlation function. Section lTTTl generalizcs 
our perturbation method to systems with many neurons 
and shows that in the asynchronous mode, only close 
neighbors contribute significantly. The final section dis- 
cusses the applications of these analytic methods to bio- 
logical systems. 



potential is the sum of all synaptic inputs. It is con- 
venient to construct our probability function in discrete 
time where, t n = nAt, n S Z and At is on the order 
of one millisecond. The probability of a spike generated 
in neuron i during time step t n is given by a threshold 
function of the membrane potential, 



Pi(t„) = Pi{Vi(t n )) 



1 



l + exp[-IM(yi(tn)-Oi)] 



(1) 



where 9i is the spike threshold and /i.; parametrizes the 
noise (/i.; ~ 1 /noise) of neuron i. The noise arises from 
ion channel fluctuations that are internal to each neu- 
ron, and from synaptic noise that arises from synaptic 
connections not explicitly included in the model. 

A presynaptic spike evoked by another neuron, j, will 
contribute to the membrane potential Vi(t n ) with a PSP 
kernel weighted by the synaptic weight, Wij. In addition 
to synaptic input, the membrane potential of neuron i is 
also given a bias (or background input), Ui(t n ), that pro- 
vides a spontaneous discharge rate or driving input to the 
neuron i. If <S* re (t m ) represents a spike train (S^ re (t m ) 
— lor at each time step) from neuron j, then the mem- 
brane potential of neuron i is, 

VS(f„) = Ui(t n ) +J2wije(t n - t m )Sy e {t m ). (2) 

The j-sum is over all neurons that have synaptic contact 
onto neuron i. To compute the ensemble average of the 
membrane potential, we use the probability functions of 
the presynaptic neurons, 



(Vi(t n )) = Ui(t n ) + ^w ij e{t n 



tm)Pj(tm) (3) 



A. Network of spike response neurons 

The model neurons used in this study are based on 
the Spike- Response Model Q . The basic idea is to con- 
struct a kernel function that describes the response of 
a biological neuron to spikes. The membrane time con- 
stant and other intrinsic properties of neurons are con- 
tained in the response function 0. In the following, 
we restrict the model to respond linearly to presynaptic 
spikes with a simple exponential decay represented by the 
postsynaptic-potential (PSP) kernel, e(t). If t pre is the 
time of the presynaptic spike, then e(t) = aexp[— a(t — 
tpre)] for t > t pre , otherwise e(t) = 0. The parameter a 
is used to describe the membrane time constant, and the 
kernel is normalized such that Je(t)dt = 1. This choice 
of kernel implies that the model is formally equivalent 
[Tfij l to integrate-and-fire models^, 0|. These kernels 
could also contain a synaptic delay, but in the examples 
we will assume that the delay is negligible. 

Spikes are generated in each model neuron i by a prob- 
ability function, Pi(t), that is dependent on a variable 
called the membrane potential, Vi(t). The membrane 



because Pj(t m ) = (Sj re (t m )). Since the spike probability 
of neuron j is functionally dependent on the spike prob- 
ability of neuron i, via its membrane potential, then in 
recurrent circuits Vi(t n ) will appear nonlinear ly on both 
sides of the equation and we must seek alternative meth- 
ods to solve for it. 



B. Fixed-point behavior and time scales 

The neurons of neural systems have been found to have 
preferred spike rates due to adaptive mechanisms that 
drive these systems to these rates and help maintain sta- 
bility of the overall system 0, [2(j • Both theoretical 
[2ll I22I |23| and experimental studies 0] have suggested 
that the establishment of a fixed point in the neural dy- 
namics is on a time scale that is much longer than the 
dynamics driven by spike responses. The presence of a 
fixed point allows one to choose a mean membrane poten- 
tial to expand the probability function in Eq. Let p 
be the stable fixed point of the long time scale dynamics: 



Pi(t r . 



P l {V i ) as t„ — > 00 



(4) 
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By expanding the probability function about this fixed 
point, we can extract a linear relation to expose the 
membrane-potential function on the right side of Eq. © 



(5) 



where Vj is the membrane potential at the fixed point. 
In the following, we will assume that the system is near 
the fixed point, p, that is independent of time. This as- 
sumption will simplify the analysis and, in one example, 
we will allow the background input, Ui(t n ), to vary in 
time. 

In the next section we will pursue this expansion by 
simplifying our network to the case with only two synap- 
tically coupled neurons to illustrate the method. How- 
ever, the technique generalizes to networks of multiple 
neurons as will be shown Section III. 



network, let Ui(t n ) be the background membrane poten- 
tial of neuron i, the membrane potential in the absence of 
synaptic connections with neuron j. Here, the symbols i 
and j are fixed labels for each of the two neurons so that 
we may easily generalize the method to larger networks 
in the following sections. The ensemble average mem- 
brane potential of neuron i is given by Eq. ©. However, 
since Pj(t m ) is dependent on the activity of neuron i due 
to the recurrent synaptic connection, we must expand 
about the background membrane potential of neuron j 
as in Eq. (j3J). The membrane potential of neuron j is 
also dependent on the spike activity of neuron i, so that 
we now have, 



II. PAIR OF COUPLED NEURONS 



Let Wij £ K be the weight (efficacy) of a synaptic con- 
nection from neuron j onto neuron i. For our 2-neuron 



J 



(Vi(t n )) = Ui(t n ) 



■ Wijpje * [1 + -Pj)(Uj(t m ) 



Vjf)] +W^ijPj{l -Pj)w. 



e (2) *P,(*„) 



(6) 



where we have introduced the notation for the multi- 
ple convolutions, e<^ * f(t n ) = ^ • • • E lK e (^-i ~ 
U K ) ■ ■ ■ t(t n — The last term in Eq. repre- 

sents the recurrent effects of neurons ?'s activity on itself 
via neuron j. Expanding the spike probability function 
of neuron i, Pi(t n ), generates another term with a factor 



J 



that contains the spike probability function of neuron j. 

For simplicity, let pi = pj — p, fii = fij = p,, 6i = 6j = 
9, and collect the terms in powers of the synaptic weight, 
Wij = Wji = w, to arrive at an expression for the average 
membrane potential, 



(Vi(t n )) = Ui(t n ) + w P J2 - P)) K ^[l + M(l - * (U a (t n ) - V a )} 

I 



(7) 



K=l 



where 



i if K is even 
j if K is odd 



(8) 



We now have an expression that is dependent on the 
background activity of either neuron, that is, we may 
calculate the effect of the recurrent connections on any 
specific neuron to an accuracy of order wp,p(l — p). 

The spike probability can also be expressed as a power 
series by substituting the expansion of the average mem- 
brane potential (Eq. 0) back into the expression for the 



spike probability function Eq. JSJ): 



p t (t n ) = £ w(i-p))^ * pV(t n ) 

K=0 



(9) 



where P^ (t n ) = P a (U a (t n )) is the background spike 
probability of neuron a, and we have defined e^(t n ) — 
S(t n ). Notice that the fixed-point, p, is still present in 
this expression. If the expansion were taken about the 
background spike probability, the expression is the same 
unless the background spike probability is not constant. 
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FIG. 2: Loop Expansion. Each term in the series expansion 
for two coupled neurons (Eq. is represented by a diagram 
of recurrent loops with K synaptic links. 



For a sensible prediction of the spike probability func- 
tion, the expansion must be bounded on the interval 
[0,1]. Because the expansion parameter, w[ip{\ — p), is 
dependent on the noise and synaptic weight, we must 
investigate the conditions for convergence of this series. 
The product, p(l — p) < 0.25 for all p, so the important 
parameters for convergence are the noise parameter and 
synaptic weight. Since \x is small for noisy networks, and 
w is small for weakly coupled networks, then the per- 
turbation series converges in the weak and noisy limit of 
these network parameters provided that * P^(t n ) is 
bounded as K — > oo. 



A. Diagrammatic methods 

The expressions for the spike probability function and 
average membrane potential can become unwieldy when 
many synaptic connections are involved with different 
weights and different neuronal types. Therefore, we in- 
troduce a book-keeping device based on Feynman-like 
diagrams that will help keep track of the terms of the 
expansion [2H I2H] . 

A chain of synaptic links is the number of synaptic 
connections that a signal passes from the starting neuron 
to the final neuron. For instance, in our example of the 
pair of coupled neurons, the second term in Fig. is 
a chain of 1 synaptic link, the third term is a chain of 
2 synaptic links. Each unique chain of synaptic links 
contributes a term to the perturbation expansion. 

For a chain of K synaptic links: 
1^2->—i->j + 1, 
the influence of the first neuron in the chain on the spike 
probability of neuron [K + 1) is given by a term contain- 
ing the following factors: 

1. WjiHjpj(l — Pj) for each synaptic link, i — > j, 

2. eW*Pf(i„), 

where 1 is the label of the first neuron in the synaptic 
chain. 

The perturbation expansion for the average membrane 
potential is handled similarly. For a chain of K synaptic 
links: 

l->2->---i->j->---K->K + l, 
the influence of the first neuron in the chain on the mem- 
brane potential of neuron (A" + 1) if given by a term 
containing the following factors: 

1. WjiPi for each synaptic link, i — > J, 



2. — pi) for each traverse, — > i — 

3. eW*[l + p(l-Pi)(tfi(*n)-Vi)]. 

The resultant perturbation series, called a loop expan- 
sion, approximates the effects of neighboring neurons by 
expanding in the order of loops. Fig. (JSJ shows the graph- 
ical representation of the first four terms of Eq. @ • The 
first term is the background spike probability, and the 
second term represents the effect of the second neuron 
on the first if there were no recurrent connections. The 
third term is then the effect of the first neuron on itself 
via its effect on the second neuron. This technique gen- 
eralizes to any number of neurons with unique synaptic 
weights, intrinsic noise and background spike probability. 
Analytic methods using another type of converging series 
have been developed previously to study population dy- 
namics of large networks |26j . 



B. Numerical comparison 

We wish to check the range of validity of our loop 
expansion of our model network by comparing the pa- 
rameter dependency of the loop expansions with a nu- 
merical simulation of a pair of coupled neurons |48j| . In 
our simulations, we set Ui(t n ) = Uj(t n ) = 9 so that the 
background spike probability, PV (t n ) = P.- (i„) = 0.5. 
With no other inputs or adaptive mechanisms in the 
model neurons, our loop-expansion fixed-point will also 
be, pi = pj = 0.5. We choose this value for our compar- 
isons because each term of the loop expansion is weighted 
by factor that is a power of Pi{l — pi); a factor that is 
maximum for f>i = 0.5. Thus, other values for P^(t n ) 
and pi would allow the series to converge faster and im- 
prove the prediction of the loop expansion. 

No time dependent inputs were present in the sim- 
ulation, thus we compared the time average and vari- 
ance of the spike probability over the sampling time to 
the prediction of the loop expansion. In the model neu- 
rons, recurrent spikes were generated by selecting pseudo- 
random numbers from a uniform distribution with a spike 
probability calculated by the spike response model (Eq. 
^1 to simulate spike generation by a Poisson point pro- 
cess. This spike generation method allows us to compare 
the loop expansion results to simulations up to the order 
of variance of a uniform distribution. The sample time 
was N = 2 x 10 6 time steps. The spike probability cal- 
culated at each time step was averaged over the sample 
time: 



P, 



1 N 



(10) 



n=l 



with the variance calculated by 



var(Pj) 



\ 



1 N 



Pi(t n )Y 
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FIG. 3: Comparison of the loop-expansion prediction with simulation for the synaptic weight dependence of the spike proba- 
bility. The time average of the spike probability function (solid line) saturates for strong synaptic weights (error bars represent 
1 standard deviation). The loop-expansion prediction (dashed lines) for if-terms of the expansion (Eq. [§J| matches the sim- 
ulation result in the presence of strong effects on the spike probability from recurrent connections. Background input is set 
to the threshold, Ui(t n ) = Uj(t n ) = 6, the noise parameter is n = 0.002, and the number of time steps for each simulation is 
N = 2 x 10 6 . 



In the figures, the variance is shown with error bars. 

Fig.|31shows the effect of synaptic strength on the the- 
oretical spike probability, where we set the noise param- 
eter to fi = 0.002. The prediction for the first 12 terms 
of the loop-expansion is indistinguishable from the sim- 
ulation for w € [—900,600]. It is important to note that 
for the excitatory synaptic weights, the lower order ex- 
pansions are better predictors of the spike probability 
for w > 500. Since the product (ip(l - p) = 1/2000, the 
loop-expansion diverges for w > |jup(l — = 2000. 

The observed deviation of the loop-expansion result 
from the simulation is in a smaller range of the synap- 
tic weights than would be predicted from the radius 
of convergence of the loop-expansion. The reason for 
this discrepancy is that there is coupling of higher order 
moments of the spike probability distribution that were 
truncated in Eq. |SJ By adopting a linear approximation 
of the spike-probability function (Eq. ^) , we have lost in- 
formation about when the higher moments of the distri- 
bution become large, as in the case w < —900. The linear 
approximation of the spike probability function also trun- 
cates information about the saturation of the spike prob- 
ability for w > 500. On the other hand, in the region of 
synaptic weight values where the loop-expansion is valid, 
the spike probability, Pi(t n ), is in the range [0.3,0.8], 
which demonstrates a significant effect of the recurrent 
connectivity. 

We also compared the dependency of the spike prob- 
ability on the noise of the system with numerical simu- 
lations for the case where the recurrent synaptic weights 
were set at w = —500. We found that the variance 
of the spike probability covers the prediction for much 
of the range of noise values. The time average of the 



spike probability function is not strongly influenced by 
the noise parameter, and the loop-expansion prediction 
for -fT-terms of the expansion (Eq. matches the simula- 
tion when the noise is high (1/fi > 300) . The deviation of 
the prediction becomes prominent where the loop expan- 
sion diverges, the noise is low. For much of the parame- 
ter range, the higher order terms of the expansion do not 
greatly improve the prediction of spike probability. This 
lack of improvement greater predictive accuracy in the 
high noise is the basis of the mean-field approximation 
[13j where only the first two terms of the loop expansion 
are used. 



C. Stimulation effects 

We can calculate the effects of surrounding neurons 
on the response of a single neuron to a current injection 
by using our expression for the membrane potential of 
the model neuron. This type of calculation would be 
useful for quantifying the functional strength of synaptic 
connections to a neuron embedded in a network when 
the neuron can only be examined with single cell record- 
ings. If the number of synaptic connections can be es- 
timated using morphological studies, and the noise and 
baseline activity measured by a single cell recording, then 
the effects of neighboring neurons on the response of the 
recorded neuron to a step stimulus can be used to fit the 
recurrent synaptic strength. 

In the simple case of two synaptically coupled neu- 
rons, we can represent a constant current injection into 
one neuron using the following background membrane 
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FIG. 4: Comparison of the loop-expansion prediction with simulation for the synaptic weight dependence of the average 
membrane-potential loop-expansion during a pulsed, intracellular simulation. The time average of the membrane potential 
(solid line) during a repeated step-change of the background input by 20% of the baseline membrane potential (error bars 
represent 1 standard deviation) . The loop-expansion prediction (dashed lines) for if-terms of the expansion (Eq. [5J matches 
the simulation for weak synaptic weights. Parameters are as in Fig. [3] except the number of time steps for each simulation is 
7V = 2xl0 5 . 



potentials (from Eq. HJ) 
Ui(t n ) = 



V + U , for t n during the stimulus, 

V otherwise, 



Uj(t n ) = V, forall*„. 



(12) 



This definition for the background potential alters the 
our loop expansion of the membrane potential (Eq. |7J). 
Since Ui(t n ) — V = U , and Uj(t n ) — V = 0, we have 



K=l 



(13) 



r 



where we have assumed that the stimulus time is long 
enough to integrate out the EPSP kernel. 

Comparisons with numerical simulation are shown in 
Fig. As in previous simulations, we calculated the 
membrane potential, Vi(t n ), and calculated spiking of the 
neuron pair using spike response models. In the simula- 
tion, one of two coupled neurons was stimulated with a 
with a step increase of its baseline membrane potential 
by 20% so that U — 0.2V. The stimulus was presented 
every M = 200 time steps and lasted for L — 20 time 
steps for a total sample time of N — 2 x 10 5 time steps. 
The membrane potential was time averaged during the 
stimulation presentations, 



Vi 



M 



N/M L 



^v(t mM+ i), 



(14) 



m+l i=l 

and the variance was calculated by 



var(Vi) = 



\ 



M 



N/M L 



mM+l 



)) 2 - 



(15) 



m+l ;=i 



The result show that the loop expansion can predict the 
response of a neuron embedded in a network with a broad 
range of synaptic strengths. In these expressions, the 
results could be generalized to the case where toy ^ Wji 
if we recover a sum over the synaptic weights and alter 
the expressions appropriately. 

In the example presented here, we assumed that the 
stimulation time was long compared to the EPSP so that 
we integrated out the temporal dependencies for simplic- 
ity. However, the change in spike rate caused by the stim- 
ulation would provide a transient correlation that persist 
for the duration of ^ K) {t n ) jHl^l]. This would suggest 
a mechanism for persistent currents leading to possible 
memory effects. 



D. Correlation functions 

We can apply our diagrammatic method directly to 
the problem of calculating correlation functions in a net- 
work of coupled neurons. For a pair of neurons, we refer 
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FIG. 5: Comparison of corelation functions with numerical simulation results. The correlation functions from the simulation 
were computed by Eg. 1191 and average over 10 trials (solid line, error bars represent 1 standard deviation). The loop-expansion 
prediction (dashed lines) for a 6-term expansion (dashed lines) matches the simulation when the synaptic strength, w, was 
within the valid range found in Fig. [3] The recurrent synaptic contact in the simulation had a 1-step time-delay as is evident 
in the joint-correlation functions, dj . The parameters are the same as those used in Fig. [J] 



to the loop-expansion diagram (Fig. [2J, along with the 
diagrammatic rules, and choose those terms that repre- 
sent the type of correlation that we wish to calculate. 
For an auto-correlation function, we use only the even 
power terms that represent the effect of a neuron on it- 
self through the other neuron in the network, and for 
a joint-correlation function, we use only the odd power 
terms the give the influence of one neuron on the other. 

1. Parameter dependency of correlation functions 

We wish to calculate the spike correlation function (2^| , 

Cy(* n ) = (((Si(t n )) - Si(t n )) X 

«S,(tn + Sn)) -Sj(t n + S n ))}. (16) 

In our example, we have let the background membrane 
potential be a constant so that (Si(t n )) = p. To obtain 
an expression for the auto-correlation function (j = i), 
we multiply the even terms of the loop expansion (Fig. 
EJ) by an overall factor of p that arises from the product 
of spike probabilities, 

oo 

Cu(t n ) = f J2(^P(l-P)) 2K C {2K) (tn). (17) 
K=l 

The terms of the loop-expansion represent the effect of 
neuron i on itself via propagation K times around the 
loop. If the background membrane potential were not 
constant, then there would have been terms in Eq. 1171 
that result from the temporal correlations of Ui(t n ). 

The joint-correlation function is found by multiplying 
the odd terms of the loop-expansion by the average spike 
probability, 

oo 

c«(*») =f E^ 1 -? 5 )) 2 ^ 1 ^" 1 ^™)- ( 18 ) 



The first term in the series is the correlation as a result of 
the direct effect of neuron- j on neuron-i, the second term 
in the series is the effect for the signal that has traveled 
once around the loop. 

2. Numerical comparison with correlation functions 

In order to test our analytic expressions for the cor- 
relation functions of two mutually connected neurons, 
we simulated a pair of spike-response model neurons and 
computed the spike correlation functions using the for- 
mula: 

1 N 

Ct; m (t n ) = - Si(t m )Sj(t m + t n ) ~ SiSj (19) 

m— 1 

where R is the number of spikes in the spike train of 
neuron-i and Si is the time-averaged spike probability, 

Si = (l/AOE^Li 5 ^™)- The sam P le time was N = 
2 x 10 5 . In the case of auto-correlation functions (i = j), 

we set Cf m (0) = H3- 

A comparison is shown in Fig.[S]for the auto- and joint- 
correlation functions of neurons coupled by two choices of 
synaptic weights. In Fig. [S] we compare the prediction of 
the loop-expansion to the average correlation function of 
10 simulations. The prediction in the case of the stronger 
weight is shown to deviate from the simulation due to 
higher order terms containing higher order convolutions 
of the EPSP function, e {K) {t n ). 

E. Comparison with the high-temperature 
expansion 

A similar method of deriving analytic expressions 
for collective variables in statistical systems was devel- 
oped that is valid under conditions of high temperatures 



[3lL 133 . 13^ . Such an expansion takes advantage of the 
noise parameter in our spike probability function (Eq.QJ. 
This function is equivalent to the distribution function of 
fcrmionic systems where the parameter /i is proportional 
to the inverse temperature of the system. Expanding Eq. 
^near \i — yields 

Pi(t n ) = \ + \»{V i {t n )-6 i )-^{V i {t n )-e i f 

+\^{v l {t n )-e t f + ••• (20) 

As in our earlier loop-expansion near a fixed point of 
Pi(t n ), we substitute in the explicit expression for Vi(t n ) 
which contains a factor representing the spike probability 
of the neurons connected to neuron-i. The spike probabil- 
ity functions for the coupled neurons are also expanded, 
and the linear terms are collected so that we arrive at an 
expression for the spike probability of neuron-i. This ex- 
pansion is the sum of its spike probability in the absence 
of synaptic connections plus a series of perturbations due 
to recurrent connections, 

Pi(tn) = if (*») + E (\^) K Z {K) * if (U (21) 
K=l 

where a = if K is even, and a = 1 if K is odd. The ex- 
pression obtained from the high-temperature expansion 
is therefore equivalent to the fixed-point loop-expansion 
for the special case when the fixed-point p — 1/2. Thus, 
the fixed-point loop-expansion can be thought of as a 
high-temperature expansion with a bias and with a tem- 
poral structure in the interaction kernels. 

Although the fixed-point loop-expansion was devel- 
oped for biological systems that possess homeostatic 
adaptation mechanisms, there are advantages to this con- 
nection with the high-temperature expansion. First, for 
more complicated networks, the diagrams that weight the 
terms of the perturbation expansion have already been 
enumerated in the literature j3^|- Second, the connec- 
tion may be made with other techniques from statistical 
mechanics |35[ to calculate properties of neural networks 
that could be used to test models of biological neural 
systems. 

III. NETWORKS WITH MULTIPLE NEURONS 

A. Effective spike propagation 
Diminishing efficacy and equivalent networks 

The propagation of spikes across multiple neurons in a 
neural network can be estimated by the loop-expansion 
of the correlation function for a periodic chain of neu- 
rons [6j, where each neuron is recurrently coupled to its 
nearest neighbors and the first neuron is coupled to the 
last. When more than two neurons are on a network, the 
number of loop diagrams that contribute to each term of 
the loop-expansion must be enumerated. 



S 




FIG. 6: Loop diagrams for the K = 5 term in the correlation 
function, Cij(t n ), of a periodic chain of neurons. There are 5 
synaptic links with two loops. The computation is symmetric 
under the transposition of i and j. 

Consider the joint-correlation function between a 
neighboring pair of neurons in the network depicted in 
Fig. n The terms of the correlation function (Eq. I18|) 
must by weighted by the number of distinct loop dia- 
grams, A l K that can be constructed with (K — l)/2 loops, 
where / is the number of synaptic links that separate the 
neurons (j — i ± I). Then we have 

oo 

C ij (t n )=p 2 J2 A l K (u^p(l-p)f K - 1 e^ K - 1 \t n ).(22) 
K=l 

For nearest neighbors (j = i + 1), the coefficients are: 
A\ = 1, Aq = 3, A\ = 8. The distinct loops contributing 
to A$ are shown in Fig. 

We compared the prediction of the loop-expansion for 
the joint-correlation function with a simulated periodic 
chain network of 10 neurons (Fig. [7J. The model pa- 
rameters were the same as used to generate Fig. |21 with 
the synaptic weights set at w = —500. We averaged 
over 10 simulations and the number of time steps in each 
simulation was N = 2 x 10 5 . The results are shown 
in Fig. where the joint-correlation function was com- 
puted in the simulation for 5 different locations. The 
loop-expansion was computed for the nearest neighbor 
(j = i + 1), and the neuron two steps away (j = i + 2). 
The joint-correlation function of the neuron pair that are 
separated by 2 synaptic links (j — i + 2), has coefficients: 
A\ = 1, A\ = 4, A\ = 12. The joint-correlation functions 
for neuron pairs separated by more than 2 synaptic links 
were indistinguishable from the noise. 

The result that neurons located further than two steps 
away had a vanishing correlation function suggests a 
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FIG. 7: Diminishing efficacy and equivalent networks. The joint-correlation function of neuron pairs that are separated by 
more than 2 synaptic connections are indistinguishable from noise. A comparison of joint-correlation functions with numerical 
simulation results for Cj^+i and d,i+2 demonstrate that the loop-expansion generalizes to networks of multiple neurons. The 
parameters are the same as those used in Fig. [3]witli the synaptic weights set at w = —500. 



rapid, diminishing efficacy of one neuron's synaptic ac- 
tion on other neurons in the network. This diminishing 
efficacy is in a parameter range where network recurrence 
reduces the spike probability of each neuron by 20% (see 
Fig. |2| . The signals that may begin with one of the 
neurons do not propagate far, even though the local con- 
nectivity has a strong effect on each individual neuron. 

Since the correlation function is indistinguishable from 
zero after 3 synaptic links, we could have derived identi- 
cal results with a periodic chain network of only 5 neu- 
rons. Thus large networks may be simulated by smaller, 
equivalent networks to yield the same results for neuron 
activity variables such as spike probability and correla- 
tion functions. 



B. Continuum limit of spatial components 

In large biological neural networks, identifying and an- 
alyzing the multiple connections between a large number 
of neurons neurons can be cumbersome. A mean-field 
approach has been developed 0, |3(| to deal with large 
populations of interacting neurons where the neural pop- 
ulations are represented by a continuous field. In our no- 
tation, the spike-probability function will be generalized 
to include spatial components, P(x, t n ), where x G M" 



J 



and 1 < n < 3. The spike probability at each point in 
the network has the same dependence on the generalized 
membrane potential, V(x, t n ), as in the discrete neuron 
case (Eq. nj, but now the neurons are labeled by their 
spatial location, x. 

Synaptic interactions are introduced by extending the 
synaptic weights to a synaptic density, tw(x), so Eq. |3| 
becomes 

(V(x,t n )) - U(x,t n ) (23) 
+ / dx'^u;(x-x')e(t ro -t n )P(x',t m ) 

J ro 

We can combine the synaptic density with the PSP ker- 
nel to define an effective PSP kernel that depends on 
both space and time, W(x,t n ) = w(x)e(t n ). Repeating 
our calculations from the discrete case (Section [nj we 
arrive at the loop expansion for the spike probability of 
a continuous neural network 

oo 

P(x,i n ) = £>P(1 -p)) K W^*P u (x,t n ) (24) 

K=0 

where the convolution is now defined over spatial as well 
as temporal variables, 



W {K) * f{x,t n ) = J dx.i ^2 ■ ■ ■ J dx K ^w{x K _ x - x K )e(t lK _ 1 -t iK )- ■■w{x - x x )e{t n ~ t il )f{x,t il ), (25) 



r 



where all n components of x are integrated. 

The loop corrections to any mean-field model of a large 
neural network can now be computed for the case where 
each location in the network is driven by a different in- 



put as represented by the background spike probability, 
P u {x,t n ). Although we have only included one type of 
neuron in the above calculations, we may further gener- 
alize the PSP kernel to couple different types of neurons, 
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as we will demonstrate in the next section. 



IV. APPLICATIONS TO BIOLOGICAL 
NETWORKS 

An important application of the techniques presented 
here is to biological neural networks found in laminar 
structures with lateral synaptic connectivity such as the 
mammalian cerebral cortex and the cerebellum. In the 
cerebral cortex, pyramidal cells are coupled with lateral 
excitatory synapses. In addition, there are inhibitory in- 
terneurons that are excited by, and inhibit, pyramidal 
cells. Interacting pools of excitatory and inhibitory neu- 
rons have been modeled using mean-field approaches [3^ | 
and have been used to study instabilities under patho- 
logical conditions in the visual cortex [3^|. We will ex- 
tend these previous results because we include the time- 
dependence in the synaptic kernel function and we add 
perturbative corrections due to recurrence. 

The synaptic densities of the two neural populations 
are modeled with a Gaussian function |38|. 

Wa (x) = ^ e -^) 2 , (26) 
where a denotes whether the synapses are excitatory (a = 



Since each convolution of Gaussian functions increases 
the standard deviation, each higher order term of the 
loop extends the reach of the synaptic densities. Thus, 
as K increases in the loop expansion, each term predicts 
the distance across the cortex that neurons contribute to 
the spike probability of a given neuron. This expression 
for the spike-probability function could be improved by 
including more details about the inhibitory interneurons, 
such as by taking care to include an independent fixed 
point spike probability, p, for each neuronal type. 

The validity of the loop expansion depends on whether 
the dynamical solution is stable in both the time and 
space domain for physiological parameter settings of the 
synaptic weights and noise parameter. Previous studies 
of cortical dynamics have suggested that spatial instabil- 
ities are the result of a pathological state of the system 
[38j . We may therefore assume that, under normal phys- 
iological conditions, the loop expansion provides a valid, 
analytic expression for the neural dynamics of the cere- 
bral cortex. 

A possible use of this calculation would be to calcu- 
late the activity of pyramidal cells in the primary visual 



E) or inhibitory (a = I), w a scales the synaptic strength, 
and a a is the lateral extent of the synaptic coupling. In 
the 2-dimensional cerebral cortex, x£l 2 and x 2 = x • x. 

The effect of excitatory pyramidal cells on other pyra- 
midal cells is both excitatory and via interneurons, so we 
may eliminate explicit reference to the inhibitory popu- 
lation by incorporating the inhibitory interneurons into 
the definition of the excitatory synaptic density, 

W(x,t B )=W B (x,t B )-W>(x,t„). (27) 

where the first term represents the excitatory component 
of the synaptic density, and the second term represents 
the inhibitory component. The time dependency of the 
inhibitory PSP kernel can be approximated by convolving 
two PSP kernels because the pyramidal cells inhibit other 
pyramidal cells through 2 synaptic links. Thus, the form 
of the synaptic density becomes 

W(x,t n ) = w E ( X )e(t n ) - Wl ( X )eW(t n ). (28) 

Using this expression for the synaptic density, we may 
compute the loop expansion for the spike-probability 
function (Eq.^JJl, 



(29) 

I 

cortex that are driven by sensory information from the 
lateral geniculate nucleus [33, represented by P u (x,t n ). 
The prediction of Eq. 1291 combines the visual input with 
the lateral connections within the primary visual cortex. 
Discrepancies between the predicted spike activity and 
spike activity measured in experiments would suggest 
the influence of feedback recurrence from other regions 
of the visual cortex Therefore, the loop-expansion 
technique could be a theoretical tool to study the con- 
sequances of synaptic communication between cortical 
regions. 

A computation could also be made to predict the spike 
correlation between neurons as a function of their sepa- 
ration. Since this analytic technique provides the pa- 
rameter dependencies of correlation functions, it is pos- 
sible to deduce the effective synaptic coupling between 
neurons using data that measures the lateral spread of 
pyramidal cell axons and pair recordings to measure the 
space dependent correlation functions. These calcula- 
tions would reveal what dynamical state of the cortex 
(coherent or asynchronous) was present under physiologi- 
cal conditions, but such calculations are beyond the scope 
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of the present study. 

Another application to biological networks is to pre- 
dict neural dynamics in the presence of recurrent feed- 
back loops, as found in the mammalian thalamic-cortical 
system and the cerebellum. In the cerebellum, the mossy 
fiber input to the system is controlled by recurrent inhibi- 
tion from inhibitory interneurons that appear to control 
the processing of information carried by mossy fibers. If 
the recurrent inhibition plays a modulatory role that can 
be modeled as a perturbation of the system, then this 
technique could be valuable to investigate the effects on 
sensory processing in the granule cell layer of the cere- 
bellum gjj. 

A class of biological neural networks that falls outside 
the validity of the loop-expansion as presented here is the 
pattern generators such as those found in invertebrate 
motor systems and the spinal chord. Since these neural 
systems operate in a bursting, coherent state, the pertur- 
bation expansion would be a poor predictor of the neural 
dynamics. Techniques from dynamical systems would be 
more appropriate for the study of these systems. 

Since the loop-expansion technique is valid when the 
recurrence plays a modulatory role in a weak and noisy 
system, the method may be generalized to analyze modu- 
latory biochemical networks. The localized chemical con- 
centrations could be represented by neurons and the ki- 
netics of the biochemistry represented by PSP kernels. 
The results could predict dynamics of complex biochem- 
ical networks where many components interact. 

In summary, we have introduced an analytic method 
to compute dynamical variables that predict the spiking 
activity of neurons embedded in recurrent neural net- 
works. The method relies on an expansion of recurrent 
synaptic connections and is valid for neural networks in 
the asynchronous state where there are no global insta- 
bilities. A diagrammatic method has been introduced to 
help construct the terms of the loop-expansion, and the 
results have been compared with simulations of spiking 
neurons. The loop-expansion technique is designed to be 
applied to biological neural networks of spiking neurons, 
and includes the temporal aspects of synaptic transmis- 
sion and spike generation. Finally, we demonstrated how 
the loop-expansion technique can be applied to neural 
dynamics in the cerebral cortex to extend previous mod- 
eling studies. 
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APPENDIX A: SPIKE PROBABILITY 
FUNCTIONS 

In order to develop a mechanistic understanding of 
noise in our threshold model of a neuron, we will com- 
pare three models on neuronal variability. The objective 
of our model presented in Section llAl is to separate the 
correlated synaptic inputs from uncorrelated noise. The 
uncorrelated noise consists of random excitatory and in- 
hibitory synaptic inputs, as well as channel noise. These 
inputs are uncorrelated with the inputs described by the 
function Vi(t n ) in Eq. Q. The model neuron described 
by Eq. is an approx imation of the Gaussian limit of 
the Stein model [4l|,|42( that is developed from a stochas- 
tic differential equation for the membrane potential v(t), 



dt 



v(t) 



v(t) 



(Al) 



where r is the membrane time constant and rep- 
resents uncorrelated Gaussian noise with zero mean: 
(£(*)) = and (£(*)£(t')> = cr 2 S(t-t'). In Stein's model, 
£(t) represents inputs from balanced excitatory and in- 
hibitory synaptic noise. 

In the diffusion limit of many small synaptic inputs, 
it can be shown |43| that Stein's model approximates an 
Ornstein-Uhlenbeck process 0] with a stationary mem- 
brane potential that takes the form of a Gaussian func- 
tion, 



a) 



1 



: exp 



V 



(A2) 



where V is the stationary mean of v(t). A neuron with a 
membrane potential described by this Gaussian distribu- 
tion, and a spike threshold of 9, has the spike probability 
function, 



P, 



Gauss 



(V,f) 



<&(v, a) dv. 



(A3) 



To compare this spike probability function with Eq. , 
we equate the slope of Pcauss(V, 0) with the slope of 
Pi{Vi) at Vi = 6 to arrive at fi = ia/y/ir. A compar- 
ison of the functional form of Eq. (|A3|I and Eq. £Q) is 
shown in Fig. [S] In our treatment, we embed the ef- 
fects of uncorrelated synaptic noise as background noise 
into the spike probability function (Eq. P) . The source of 
noise parametrized by fi is uncorrelated with the synap- 
tic input from other neurons explicitly represented in the 
network and with the background input, (£(t)Ui(t')) = 0. 
Improvements to this approximation could include more 
details of how s ynap tic dynamics affect the characteris- 
tics of the noise |45j . 

It is instructive to compare the spike probability func- 
tion developed by Gerstner 0,0, because Gerstner's 
model has an explicit expression for the instantaneous 
spike-rate function, p(V), based on an Arrhenius hazard 
function |46j . Our spike probability function given by 
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FIG. 8: Comparison of three forms of spike probability functions. The solid trace is the function defined in Eq. the dotted 
trace is the function defined in Eq. 1A3L the dashed trace is the function defined in Eq. 1A4L The noise parameters a = 0.5, 
p = 4<r/y / 7r, and (3 — p/ ln(4), and the threshold, 8 — 0. The differences between these models are within the variability of the 
simulations compared with the loop expansion in previous figures. 



Eq. estimates the probability of a spike in the in- 
terval At. The sigmoid threshold function is similar to 
Gerstner's interval spike probability function found by 
integrating the probability of surviving without a spike, 
cxp(— tp(V)), over a time interval. The instantaneous 
spike-rate is p(V) = roexp[/3(V — 8)], where (3 is the 
noise parameter and rrj is the spike rate at threshold. 
Integrating exp(— tp(V)) over the interval At yields the 
spike probability function 01 

P G erstner(V,0) = 1 - CXp [-r At exp[/3(F - 9)}} . (A4) 



Equating this function with our spike probability func- 
tions at V = (assuming At = 1 ms) allows us to 
compute tq = ln(2). If we expand Pcerstner(V,9) about 
V = 6, we find that the relationship between the noise 
functions is p = ln(4)/3. Gerstner's spike probability 
function is compared with the other two probability func- 
tions in Fig. [S] 
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